In [2]:
%matplotlib inline
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
In [3]:
protein = pd.read_csv("nuclear.csv")
In [4]:
# mean & stdev of each class
print('Mean(cm): \n', protein.groupby(['class']).mean())
print('-'*80,'\n')
print('Stdev(cm): \n', protein.groupby(['class']).std())
print('-'*80)
In [5]:
proteinData = [protein.iloc[:,1:18], protein.iloc[:,11:21], protein.iloc[:,21:31], protein.iloc[:,31:41],
protein.iloc[:,41:51], protein.iloc[:,51:61], protein.iloc[:,61:71], protein.iloc[:,71:78]]
#print(proteinData)
In [6]:
proteinClass = protein.iloc[:, 78:82]
#print(proteinClass)
In [7]:
# 挑出要分析的 attribute，和 class
sub_protein = pd.concat([proteinData[0],proteinClass], axis=1)
sub_protein.describe(percentiles=[])
Out[7]:
In [8]:
na_cols = sub_protein.columns[sub_protein.isna().any()].tolist()
print(na_cols)
In [9]:
## drop
proteinClean = sub_protein.dropna()
#print(proteinClean)
In [10]:
## fill missing value with mean of each group
fill_protein = sub_protein.copy()
for n in na_cols:
fill_protein[n] = fill_protein.groupby(['class'], sort=False)[n].apply(lambda x: x.fillna(x.mean()))
#protein[x] = protein.groupby("class").transform(lambda x: x.fillna(x.mean()))
#for c in proteinClass['class'].unique():
#print(c)
#proteinFill = protein.fillna(,inplace=True)
In [11]:
# 印出有 missing value 的 column
## 未處理 missing value 的原資料
na_cols = sub_protein.columns[sub_protein.isna().any()].tolist()
print(na_cols)
## fill_protein 已用平均值填補完，沒有 missing value，印出 empty list
na_cols = fill_protein.columns[fill_protein.isna().any()].tolist()
print(na_cols)
In [12]:
sub_protein.describe(percentiles=[], include='all')
Out[12]:
1. 線圖¶
- 從 77 個蛋白質表現隨意挑一些出來畫折線圖
- 此例我挑前五個 attribute，['DYRK1A_N', 'ITSN1_N', 'BDNF_N', 'NR1_N', 'NR2A_N']
- 可知每個 attribute 值域分佈的範圍不同，此 dataset 沒有對他們做標準化，若要做後續分析可能要考慮這件事。
- y 軸是其數值（表現量）
- x 軸是第幾筆資料 (原資料有 1080 筆)
In [12]:
print(list(fill_protein.iloc[:,0:5]))
fill_protein.iloc[:,0:5].plot()
Out[12]:
2. 盒鬚圖¶
- 以 DYRK1A_N 為例
- 將資料分組畫 boxplot 可看出各 class 的資料離散程度
- 例如 c-CS-s 就有一些離群值，如同作業三最後提到的。
In [13]:
fig, ax = plt.subplots(figsize=(8, 6))
plt.suptitle('')
fill_protein.boxplot(column=['DYRK1A_N'], by='class', ax=ax)
Out[13]:
3. 直方圖¶
- 此例我挑前 4 個 attribute，['DYRK1A_N', 'ITSN1_N', 'BDNF_N', 'NR1_N']
- 和 box plot 一樣可看出資料離散程度，像 DYRK1A_N 和 ITSN1_N 都有一些右邊的離群值
- histogram 看不出 box plot 能表示的四分位距和平均值
- 但可看出眾數大概的位置，如 DYRK1A_N 的眾數大概落在 0.25
- (box plot 看不出來)
- y 軸是資料筆數
- x 軸是蛋白質表現程度的數值 (每個區間寬度是 0.25)
In [14]:
fill_protein.iloc[:,0:4].hist(figsize=(8,6))
Out[14]:
- 只觀察 DYRK1A_N， 將其以 Genotype 分為 Control(正常) 和 Ts65Dn(唐氏症) 分別作圖
- 可以看到 Control 有較多導致右偏的離群值，應為 c-CS-s 類的那幾個點。
In [15]:
fill_protein.groupby('Genotype').hist(column=['DYRK1A_N'])
Out[15]:
4. 散佈圖¶
- 將前 17 個 attribute 任兩個取起來畫圖 (pairplot)
- 可以看出其相關性
- 對角線則是該 attribute 分組的折線圖（類似 histogram）
In [16]:
sns.pairplot(fill_protein, hue='class')
Out[16]:
- 各 attribute 的 histogram (參考)
In [17]:
fill_protein.plot(kind='hist')
Out[17]:
例 4-1¶
- 挑出 NR1_N 和 pNR1_N 畫散佈圖，可觀察知這兩種蛋白質的表現程度成高度正相關(對於所有樣本)。
- 實際上，兩種是類似的蛋白質(皆為 NMDA 受器)，因此生物機制上的確有高度相關性。
In [13]:
Relation = pd.concat([fill_protein['NR1_N'], fill_protein['pNR1_N']], axis=1)
Relation = pd.concat([Relation, proteinClass], axis=1)
sns.pairplot(Relation, hue='class')
Out[13]:
例 4-2¶
- 一次畫上 8 種 class 的點有點繁雜
附上以基因型(Genotype)、藥物注射(Treatment)、學習刺激(Behavior)分別著色的散佈圖。
In [19]:
sns.pairplot(Relation, hue='Genotype')
sns.pairplot(Relation, hue='Treatment', palette="husl")
sns.pairplot(Relation, hue='Behavior', palette="Set2")
Out[19]:
例 4-3¶
- 挑出 DYRK1A_N 和 ITSN1_N 畫散佈圖，除了兩個 attribute 成正相關。
- 同時可發現綠色的點 (c-CS-s) 對於此兩 attribute 都算是離群值。
In [20]:
Relation2 = pd.concat([fill_protein['DYRK1A_N'], fill_protein['ITSN1_N']], axis=1)
Relation2 = pd.concat([Relation2, proteinClass], axis=1)
sns.pairplot(Relation2, hue='class')
Out[20]:
5. 不使用圓餅圖¶
- 以我個人感受而言，我覺得圓餅圖不適合作為資料視覺化
人眼對於二維圖形的大小比例關係並無法判斷的很精準， 也許判斷相對大小沒問題，但判斷絕對大小卻很困難。
google "pie chart" 其實也可看到很多爭論。
- 參考：https://blog.funnel.io/why-we-dont-use-pie-charts-and-some-tips-on-better-data-visualizations
1. IQR¶
- 以前 17 個蛋白質表現量為例，其實直接下 describe() 就會列出四分位數 (25%, 50%, 75%)
In [21]:
fill_protein.describe()
Out[21]:
- 或我們可以只取 25th 百分位數和 75th 百分位數求差，即為四分位距 (IQR)
In [22]:
Q1 = fill_protein.quantile(q=0.25)
Q3 = fill_protein.quantile(q=0.75)
IQR = Q3-Q1
row_name=['Q3', 'Q1', 'IQR']
IQR_DF = pd.DataFrame([Q3, Q1, IQR], row_name)
IQR_DF.assign()
Out[22]:
2. Sorting¶
- 取前五個蛋白質表現量為例，預設為 descending
- 為了方便看結果，下 head() 只取前五 row
In [23]:
fill_protein.iloc[:,0:5].sort_values(by=['NR1_N']).head()
Out[23]:
- 其實也可以用多變數排序，例如先排 DYRK1A_N，有相等的 DYRK1A_N 再以 ITSN1_N 排序
- 但因為此資料集為浮點數，幾乎沒什麼相同的值，所以通常用一個變數就能排序完了
In [24]:
fill_protein.iloc[:,0:5].sort_values(by=['DYRK1A_N', 'ITSN1_N']).head()
Out[24]:
3. correlation matrix¶
- 為方便看結果，挑前十項蛋白質表現量印出 correlation matrix 及做圖 (heatmap)
3-1 correlation matrix¶
In [57]:
pcorr=fill_protein.iloc[:,0:10].corr()
pcorr.assign()
Out[57]:
3-2 correlation matrix with heatmap¶
In [63]:
pcorr.style.background_gradient().set_precision(2)
Out[63]:
3-3 heatmap to show correlation¶
In [67]:
f, ax = plt.subplots(figsize=(10, 8))
sns.heatmap(pcorr, mask=np.zeros_like(pcorr, dtype=np.bool), cmap="Blues",
square=True, ax=ax)
Out[67]:
4. box plot¶
- 以 DYRK1A_N 為例
- 將資料分組畫 boxplot 可看出各 class 的資料離散程度
- 例如 c-CS-s 就有一些離群值，如同作業三和作業四提到的。
In [25]:
fig, ax = plt.subplots(figsize=(8, 6))
plt.suptitle('')
fill_protein.boxplot(column=['DYRK1A_N'], by='class', ax=ax)
Out[25]:
- 也可以用 Genotype 分別畫 box plot
In [26]:
fig, ax = plt.subplots(figsize=(8, 6))
plt.suptitle('')
fill_protein.boxplot(column=['DYRK1A_N'], by='Genotype', ax=ax)
Out[26]:
5. histogram with density¶
- 以前 4 個蛋白質表現量為例
In [99]:
sns.set(style="white", palette="muted", color_codes=True)
# Set up the matplotlib figure
f, axes = plt.subplots(2, 2, figsize=(9,9))
sns.distplot(fill_protein['DYRK1A_N'], color="b", ax=axes[0, 0])
sns.distplot(fill_protein['ITSN1_N'], color="r", ax=axes[0, 1])
sns.distplot(fill_protein['BDNF_N'], color="g", ax=axes[1, 0])
sns.distplot(fill_protein['NR1_N'], color="m", ax=axes[1, 1])
#plt.setp(axes)
#plt.tight_layout()
plt.show()
6. others discussion¶
- 這份資料的數值可解釋性我覺得不算高，不太好理解每欄的數值到底代表什麼意思，但由這些作圖其實可以找出一些分布相似或是相關性高的屬性。
- 例如 [DYRK1A_N, ITSN1_N] 或是 [NR1_N, pNR1_N]
- 這個資料集還蠻適合做 clustering 或 classification，接下來我應該會試著用 self-organizing map 來做 clustering
找出相關性高的 attribute 除了更方便分群合併相鄰 attribute 外，
也適合在降維(PCA,LDA)的時候把不相關的幾個 attribute 先剔除，降低整體複雜度。
1. 解釋為何使用這些變數當作 x 變量或 y 變量¶
- 挑選 X = NR1_N 和 y = pNR1_N，因已知兩種蛋白質為生物上同源且化學結構類似，應有高度相關性，方便迴歸預測
- 由 Hw5 使用過的 Relation DataFrame 取出這兩行資料。
- 將 X, y 取 20% 出來做 testing data，其餘 80% 為 training data
In [245]:
Relation.describe(percentiles=[])
Out[245]:
In [396]:
X, y = Relation.iloc[:,0].values.reshape(-1,1), Relation.iloc[:,1].values.reshape(-1,1)
In [397]:
from sklearn.metrics import mean_squared_error
row = ['linear', 'knn', 'ridge', 'linear^3', 'linear^4+intercept']
col = ['MSE', 'Cor']
regResult= pd.DataFrame(index=row, columns=col)
In [398]:
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)
In [423]:
from sklearn.linear_model import LinearRegression
# Train linear model by training set
reg = LinearRegression().fit(X_train, y_train)
# Make predictions using the testing set
y_pred_lm = reg.predict(X_test)
# The coefficients
print('Coefficients: \n', reg.coef_)
print('Intercept: \n', reg.intercept_)
print('linear model Correlation: \n', reg.score(X_train, y_train))
# Plot outputs
plt.scatter(X_test, y_test, color='black', label='test data')
plt.plot(X_test, y_pred_lm, color='blue', linewidth=3, label='linear model prediction')
regResult.iloc[0, 0] = mean_squared_error(y_test, y_pred_lm)
regResult.iloc[0, 1] = model1.score(X_train, y_train)
plt.legend()
plt.show()
2-2 k-nearest-neighbor(kNN) Regression with arithmetic progression¶
- 用綠線表 kNN model 預測 testing data 的值 ( $\hat{y}_{kNN}$ )
- k 取 5，其實因為這兩個變數相關性過高，導致 k 取多少其實差不多
- 有用 uniform 和 distance 兩種方式來計算 neighbor，但結果其實差不多。
- 和 linear model 做比較
- 用藍線表 linear model 預測 testing data 的值 ( $\hat{y}_{lm}$ )
- 用黑點表 testing data 的散佈圖 ( $y$ )a
並印出 kNN model 的 correlation
可以看到右上角有一個比較遠的點有點像離群值，當 test data 挑到那個點的時候會導致 knn model 有點歪掉
In [400]:
from sklearn import neighbors
# Fit regression model
k = 5
T = np.linspace(X_test.min(), X_test.max(), 108)[:, np.newaxis]
plt.figure(figsize=(12,12))
for i, weights in enumerate(['uniform', 'distance']):
knn = neighbors.KNeighborsRegressor(k, weights=weights)
y_pred_knn = knn.fit(X_train, y_train).predict(T)
plt.subplot(2, 1, i + 1)
plt.scatter(X_test, y_test, c='k', label='test data')
plt.plot(T, y_pred_knn, c='g', label='knn prediction', linewidth=4)
plt.plot(X_test, y_pred_lm, color='blue', linewidth=4, label='linear model prediction')
#plt.axis('tight')
plt.ylim(y_test.min()-y_test.min()/5, y_test.max()+y_test.max()/5)
plt.legend()
plt.title("KNeighborsRegressor (k = %i, weights = '%s')" % (k,
weights))
# The coefficients
print('kNN Correlation: \n', knn.score(X_train, y_train))
plt.tight_layout()
plt.show()
2-3 k-nearest-neighbor(kNN) Regression with testing data¶
- 用綠線表 kNN model 預測 testing data 的值 ( $\hat{y}_{kNN}$ )
- k 取 5，其實因為這兩個變數相關性過高，導致 k 取多少其實差不多
- 有用 uniform 和 distance 兩種方式來計算 neighbor，但結果其實差不多。
- 和 linear model 做比較
- 用藍線表 linear model 預測 testing data 的值 ( $\hat{y}_{lm}$ )
- 用黑點表 testing data 的散佈圖 ( $y$ )a
- 並印出 kNN model 的 correlation
- testing data 若沒有排序的話，kNN 的折線圖容易橫向的去連接下一個點，看起來不像迴歸線
In [401]:
X_train2, X_test2, y_train2, y_test2 = train_test_split(X, y, test_size=0.2)
In [402]:
from sklearn import neighbors
# Fit regression model
k = 5
T = np.linspace(X_test2.min(), X_test2.max(), 216)[:, np.newaxis]
plt.figure(figsize=(12,12))
for i, weights in enumerate(['uniform', 'distance']):
knn = neighbors.KNeighborsRegressor(k, weights=weights)
y_pred_knn = knn.fit(X_train2, y_train2).predict(X_test2)
plt.subplot(2, 1, i + 1)
plt.scatter(X_test2, y_test2, c='k', label='data')
plt.plot(X_test2, y_pred_knn, c='g', label='prediction', linewidth=3)
plt.plot(X_test, y_pred_lm, color='blue', linewidth=3)
#plt.axis('tight')
plt.ylim(y_test2.min()-y_test2.min()/5, y_test2.max()+y_test2.max()/5)
plt.legend()
plt.title("KNeighborsRegressor (k = %i, weights = '%s')" % (k,
weights))
# The coefficients
print('Correlation: \n', knn.score(X_train, y_train))
plt.tight_layout()
plt.show()
- 將 X_test 排序完之後再放入 model 去 predict，畫出來就比較好看了
In [410]:
from sklearn import neighbors
# Fit regression model
k = 5
T = np.linspace(X_test2.min(), X_test2.max(), 108)[:, np.newaxis]
plt.figure(figsize=(12,12))
for i, weights in enumerate(['uniform', 'distance']):
knn = neighbors.KNeighborsRegressor(k, weights=weights)
y_pred_knn = knn.fit(X_train2, y_train2).predict(np.sort(X_test2.reshape(1,-1)).reshape(-1,1))
plt.subplot(2, 1, i + 1)
plt.scatter(X_test2, y_test2, c='k', label='data')
plt.plot(np.sort(X_test2.reshape(1,-1)).reshape(-1,1), y_pred_knn, c='g', label='prediction', linewidth=3)
plt.plot(X_test, y_pred_lm, color='blue', linewidth=3)
#plt.axis('tight')
plt.ylim(y_test2.min()-y_test2.min()/5, y_test2.max()+y_test2.max()/5)
plt.legend()
plt.title("KNeighborsRegressor (k = %i, weights = '%s')" % (k,
weights))
regResult.iloc[1, 0] = mean_squared_error(y_test2, y_pred_knn)
regResult.iloc[1, 1] = model1.score(X_train2, y_train2)
# The coefficients
print('Correlation: \n', knn.score(X_train, y_train))
plt.tight_layout()
plt.show()
2-4 Polynomial Ridge Regression¶
- Ridge Regression
- 用預設的 Ridge regression 做 model，alpha 值預設為 1
- 加入平方項和三次方項建 model
- 紅線表示 $\hat{y}_{1} = w_1X$
- 橘線表示 $\hat{y}_{2} = w_1X+w_2X^2$
- 黃線表示 $\hat{y}_{3} = w_1X+w_2X^2+w_3X^3$
- 因變相相關性過大，其實結果跟 linear model 看起來很像，不太需要用 ridge regression
- 隨著使用的變數變多，Correlation 有隨著提升
In [412]:
from sklearn.linear_model import Ridge
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import make_pipeline
colors = ['red', 'orange', 'teal', 'yellowgreen', 'gold']
lw = 2
#plt.plot(X_train, y_train, color='cornflowerblue', linewidth=lw,
# label="ground truth")
plt.scatter(X_test2, y_test2, color='navy', s=30, marker='o', label="training points")
for count, degree in enumerate([1, 2, 3]):
model = make_pipeline(PolynomialFeatures(degree), Ridge())
model.fit(X_train2, y_train2)
y_polyRidge = model.predict(np.sort(X_test2.reshape(1,-1)).reshape(-1,1))
print('Correlation%d: \n' % degree, model.score(X_train2, y_train2))
plt.plot(np.sort(X_test2.reshape(1,-1)).reshape(-1,1), y_polyRidge, color=colors[count], linewidth=lw,
label="degree %d" % (degree))
regResult.iloc[2, 0] = mean_squared_error(y_test2, y_polyRidge)
regResult.iloc[2, 1] = model1.score(X_train2, y_train2)
plt.legend()
plt.show()
2-5 Polynomial Linear Regression¶
- 用深綠粗線表 model1 ( $\hat{y}_{m1} = w_1X + w_2X^2 + w_3X^3$ )
用淺綠細線表 model2 ( $\hat{y}_{m2} = w_0+ w_1X + w_2X^2 + w_3X^3 + w_4X^4$ )
雖然可以決定 model 要不要考慮 intercept 項，但因這兩個變數相關性太高，導致不管有沒有考慮 train 出來的 model 都長差不多。
- 因此 model2 除了增加 intercept 項，再多增加一個 4 次方項以和 model1 比較。
- 其實會發現 Correlation 並沒差多少，model 和 data 已經相當 fit 了。
In [421]:
from sklearn.preprocessing import PolynomialFeatures
from sklearn.linear_model import LinearRegression
from sklearn.pipeline import Pipeline
import numpy as np
model1 = Pipeline([('poly', PolynomialFeatures(degree=3)),
('linear', LinearRegression(fit_intercept=False))])
# fit to an order-3 polynomial data
model1 = model1.fit(X_train2, y_train2)
model1.named_steps['linear'].coef_
model1.named_steps['linear'].intercept_
# The coefficients
print('Coefficients: \n', model1.named_steps['linear'].coef_)
print('Intercept: \n', model1.named_steps['linear'].intercept_)
print('Correlation: \n', model1.score(X_train2, y_train2))
y1 = model1.predict(np.sort(X_test2.reshape(1,-1)).reshape(-1,1))
regResult.iloc[3, 0] = mean_squared_error(y_test2, y1)
regResult.iloc[3, 1] = model1.score(X_train2, y_train2)
In [417]:
model2 = Pipeline([('poly', PolynomialFeatures(degree=4)),
('linear', LinearRegression(fit_intercept=True))])
# fit to an order-3 polynomial data
model2 = model2.fit(X_train2, y_train2)
model2.named_steps['linear'].coef_
model2.named_steps['linear'].intercept_
# The coefficients
print('Coefficients: \n', model2.named_steps['linear'].coef_)
print('Intercept: \n', model2.named_steps['linear'].intercept_)
print('Correlation: \n', model2.score(X_train2, y_train2))
y2 = model2.predict(np.sort(X_test2.reshape(1,-1)).reshape(-1,1))
regResult.iloc[4, 0] = mean_squared_error(y_test2, y2)
regResult.iloc[4, 1] = model2.score(X_train2, y_train2)
In [419]:
plt.figure(figsize=(8, 6))
plt.scatter(X_test2, y_test2, color='navy', s=30, marker='o', label="training points")
plt.plot(np.sort(X_test2.reshape(1,-1)).reshape(-1,1), y1, c='teal', label='fit_intercept=True', linewidth=5)
plt.plot(np.sort(X_test2.reshape(1,-1)).reshape(-1,1), y2, c='yellowgreen', label='fit_intercept=False')
plt.legend()
plt.show()
3. Model Evaluation¶
- 使用基本的 Mean Square Error
- 結果最基本的 linear model 就已經是 error 最小的了...
In [426]:
regResult.assign()
Out[426]:
4. Correlation between X & y¶
- 可以看到這次挑的變數相關性過大，導致各種 model 結果其實差不多。
- 這份 dataset 好像沒有比較適合迴歸的資料，也許有機會再用其他筆資料試試看。
In [424]:
regResult.assign()
Out[424]: